Set up

library(assertthat)
library(reshape2)
library(ggplot2)
library(GGally)
library(plyr)
library(rootSolve)
library(deSolve)
library(knitr)
library(pander)
library(cowplot)

#sourceFiles <- 'R/invasion.R'
#l_ply(sourceFiles, source)

set.seed(100) #reproducible random selection

sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X Yosemite 10.10.5
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_0.7.0  pander_0.6.0   knitr_1.15.1   deSolve_1.14  
##  [5] rootSolve_1.7  plyr_1.8.4     GGally_1.3.0   ggplot2_2.2.0 
##  [9] reshape2_1.4.2 assertthat_0.1
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.8        magrittr_1.5       munsell_0.4.3     
##  [4] colorspace_1.3-1   stringr_1.1.0      tools_3.3.2       
##  [7] grid_3.3.2         gtable_0.2.0       htmltools_0.3.5   
## [10] yaml_2.1.14        lazyeval_0.2.0     rprojroot_1.1     
## [13] digest_0.6.10      tibble_1.2         RColorBrewer_1.1-2
## [16] evaluate_0.10      rmarkdown_1.2      stringi_1.1.2     
## [19] scales_0.4.1       backports_1.0.4    reshape_0.8.6

Basic model description

The goal of this model to simulation biologically mediated decomposition. There are two main pools, Substrate (\(C\)) and Biomass (\(B\)). Carbon is added to the substrate pool via inputs (\(I\)), leached from both pools at a rate proportional to the carbon in the pool (\(h\) fraction of the biomass pool \(B\), and \(m\) fraction of the substate pool \(C\)), and transfer between the substrate \(C\) pool and biomass \(B\) pool via a Monod uptake with a maximum rate \(v_{max}\), half-saturation constant \(k\), and carbon use effiency \(\epsilon_{cue}\).

\(\frac{dB}{dt}=\frac{\epsilon_{cue}vBC}{k+C}-hB\)

\(\frac{dC}{dt}=I-mC-\frac{vBC}{k+C}\)

At steady state (\(\frac{dB}{dt}=\frac{dC}{dt}=0\)) this implies that:

\(B=\frac{\epsilon h}{I-mC}\)

\(C=\frac{hk}{\epsilon v-h}\)

Calculate the leaching and turnover rate

Assume that the input rate, biomass, and carbon pools are given targets, and that we are exploring cue, v, and k. We can thus define the leaching and turnover rates at steady state:

find_h <- function(C, cue, v, k){
  return(C*cue*v/(k+C))
}

find_m <- function(cue, v, h, k, I, B){
  return((cue*v-h)/(cue*k)*(cue/h*I-B))
}

If we further assume that there is a trade off between both 1) uptake \(v\) and carbon use effiency (\(\epsilon\)) such that \(v=\frac{e^{b\epsilon_cue}-e^b}{1-e^b}\) and the half saturation constant (\(k\)) such that \(k=a v +k_{min} = \frac{a}{1-e^b}(e^{b\epsilon_cue}-e^b)\). Defined as:

Visualize trade-off between cue vs uptake and k

cue_v_tradeoff <- function(b, vmax, cue){
  return(vmax*(exp(b*cue)-exp(b))/(exp(0)-exp(b)))
}

tradeoff.df <- adply(.data=c(-0.1*1:9, -1*1:10, 0.1*1:9, 1:10), .margins=c(1), .id=c('id'), .fun=function(b){
  ans <- data.frame(b=b, cue=seq(0, 1, length=1000))
  ans$v <- cue_v_tradeoff(b=b, vmax=2, cue=ans$cue)
  return(ans)
})

ggplot(tradeoff.df) + geom_line(aes(x=cue, y=v, group=b, color=b))

v_k_tradeoff <- function(kmin, k_v_slope, v){
  return(kmin+k_v_slope*v)
}

tradeoff.df$k <- v_k_tradeoff(kmin=0, k_v_slope=100/2, tradeoff.df$v) #slope is max SOC over vmax range

This leads to the following steady state calculations for one biomass pool:

steadyState <- function(b, vmax, kmin, k_v_slope, h, m, I, cue){
  v <- cue_v_tradeoff(b=b, vmax=vmax, cue)
  k <- v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=v)
  
  C <- h*k/(cue*v-h)
  B <- cue/h*(I-m*C)
  return(list(B=B, C=C))
}

Steady state sensitivity to inputs

How do inputs affect stead state?

ssParm.ls <- list(b=c(-10, -1, 1, 10),
                vmax=c(0.1, 1, 10, 100),
                kmin=c(0, 1, 10, 100, 1000),
                k_v_slope=c(1/10, 1, 10, 100, 1000),
                h=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
                m=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
                cue=c(0.2, 0.4, 0.6, 0.8),
                I=c(0.1, 1, 10))


steadyState.df <- expand.grid(ssParm.ls)

#Calculate the steady state values for each parameter combination
test <- with(steadyState.df, as.data.frame(steadyState(b, vmax, kmin, k_v_slope, 
                                                       h, m, I, cue)))
steadyState.df <- cbind(steadyState.df, test)

#Check to see if it's a reasonable value over any inputs
steadyState.df$OnTarget <- with(steadyState.df, B+C > 10 & B+C < 500 & 
                                  B/(B+C)> 0.001 & B/(B+C) < 0.15)
goodPars <- ddply(steadyState.df, c('b', 'vmax', 'kmin', 'k_v_slope', 'h', 'm', 'cue'),
                  summarize, validSet=any(OnTarget))
goodPars <- subset(goodPars, validSet, -validSet)
goodPars$index <- 1:nrow(goodPars)

#Calculate the change in the steady state over different input values
ssGrad.df <- ddply(goodPars, c('b', 'vmax', 'kmin', 'k_v_slope', 'h', 'm', 'cue'),
                   function(xx){
  xx <- data.frame(xx, I=c(0.1, 1, 10)*c(rep(1, 3), rnorm(3, mean=1, sd=0.05)))
  temp <- cbind(xx, as.data.frame(with(xx, steadyState(b, vmax, kmin, k_v_slope, h, m, 
                                   I, cue))))
  ans <- temp[1:3,]
  ans$dB.dI <- (temp$B[1:3]-temp$B[4:6])/(temp$I[1:3]-temp$I[4:6])
  ans$dC.dI <- (temp$C[1:3]-temp$C[4:6])/(temp$I[1:3]-temp$I[4:6])
  return(ans)
})

ssGrad.df$OnTarget <- with(ssGrad.df, B+C > 10 & B+C < 500 & B/(B+C)> 0.001 & B/(B+C) < 0.15)

plot.df <- melt(subset(ssGrad.df, index %in% sample(unique(ssGrad.df$index), size=10),
                       select=c('index', 'I', 'B', 'C')), 
                id.vars=c('I', 'index'))
ggplot(plot.df) + geom_line(aes(x=I, y=value, color=as.factor(index))) +
  facet_wrap(~variable) +
  labs(title='Steady state carbon pool by inputs', x='Inputs', y='Pool size') +
  guides(color=FALSE)

plot.df <- melt(subset(ssGrad.df, OnTarget), 
                id.vars=c('B', 'C', 'dB.dI', 'dC.dI', 'OnTarget', 'index'))

ggplot(plot.df) + geom_boxplot(aes(x=as.factor(value), y=dB.dI)) +
  facet_wrap(~variable, scales='free') + 
  labs(title='dB/dI', x='Parameter value', y='Gradient (dB/dI)')

ggplot(plot.df) + geom_boxplot(aes(x=as.factor(value), y=dC.dI)) +
  facet_wrap(~variable, scales='free') + 
  labs(title='dC/dI', x='Parmaeter value', y='Gradient (dC/dI)')

Optimum carbon use effiency

We want the optimum carbon use effiency where biomass is maximum as steady state but that is not possible to solve analytically (you end up with \(y=xe^x\), solve for x).

optimum_cue <- function(b, vmax, kmin, k_v_slope, h, m, I){
  cue <- seq(0, 1, length=100)
  ans <- steadyState(b, vmax, kmin, k_v_slope, h, m, I, cue)
  ans$B[ans$B <= 0 | ans$C <= 0] <- -Inf
  
  best <- which.max(ans$B)
  return(list(B=ans$B[best], C=ans$C[best], cue=cue[best]))
}

Generate reasonable parameter sets without compitition

Given that we have input targets between 0.1 and 10 mg-C per g-soil per day, total carbon pools of between 10 to 500 mg-C per g-soil, and a biomass to total carbon ratio of between 0.1 and 15 percent. Let’s generate some reasonable parameters to draw on.

parm.ls <- list(b=c(-10, -1, -0.1, 0.1, 1, 10),
                vmax=c(0.1, 0.5, 1, 5, 10, 100),
                kmin=c(0, 1, 10, 100, 1000),
                k_v_slope=c(1/10, 1, 10, 100, 1000),
                h=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
                m=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
                I=c(0.1, 1, 10))
                
optParm.df <- expand.grid(parm.ls)

optParm.df <- ddply(optParm.df, names(optParm.df), function(xx){
  ans <- as.data.frame(optimum_cue(b=xx$b, vmax=xx$vmax, kmin=xx$kmin, k_v_slope=xx$k_v_slope, 
              h=xx$h, m=xx$m, I=xx$I))
  return(ans)
})

optParm.df$OnTarget <- with(optParm.df, B+C > 10 & B+C < 500 & B/(B+C)> 0.001 & B/(B+C) < 0.15)

#parm.allI <- ddply(subset(optParm.df, OnTarget), setdiff(names(parm.ls), 'I'), summarize, 
#                 numI=length(I), minI=min(I), maxI=max(I))
#ggpairs(subset(optParm.df[,c(names(parm.ls), 'OnTarget')]) 

ggplot(subset(optParm.df, OnTarget)) + geom_histogram(aes(x=cue))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(subset(optParm.df, OnTarget)) + geom_violin(aes(x=as.factor(I), y=cue)) 

#ggplot(subset(optParm.df, OnTarget)) + geom_jitter(aes(x=I, y=B+C)) + scale_x_log10()

Competition model co-existance

Consider the senerio where the native popuation is invated by a second competiting population that has a different carbon use effiency but is otherwise equivalent. Let the strategic carbon use effiency be where the native popuation is never invadable no matter what carbon use effiency value the invader has.

In otherwords both popuations are stable where one is either 0 or \(\frac{\epsilon_n v_n}{k_n + C} = \frac{\epsilon_i v_i}{k_i + C}\) where \(n\) and \(i\) stand for the native and invading populations respectively.

Again, not solvable analytically given our trade-off function but numerically we can calculate the following:

strategic_cueI <- function(b, vmax, kmin, k_v_slope, m, h, I){
  cue <- seq(0, 1, length=100)
  v <- cue_v_tradeoff(b=b, vmax=vmax, cue)
  k <- v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=v)
  
  ans <- steadyState(b, vmax, kmin, k_v_slope, h, m, I, cue)
  
  measure <- cue*v/(k+ans$C)
  measure[measure <= 0] <- NA
  flag <- which.max(measure)
  
  return(list(cue_strat=cue[flag], C_strat=ans$C[flag], B_strat=ans$B[flag]))
}

stratParm.df <- optParm.df
stratParm.df <- rename(stratParm.df, c('cue'='cue_opt', 'B'='B_opt', 'C'='C_opt', 'OnTarget'='validPools_opt'))
stratParm.df <- ddply(stratParm.df, names(stratParm.df), function(xx){
  return(as.data.frame(strategic_cueI(b=xx$b, vmax=xx$vmax, kmin=xx$kmin, 
                        k_v_slope=xx$k_v_slope, m=xx$m, h=xx$h, I=xx$I)))
})

stratParm.df$validPools_strat <- with(stratParm.df, B_strat+C_strat > 10 & 
                                        B_strat+C_strat < 500 & 
                                        B_strat/(B_strat+C_strat)> 0.001 & 
                                        B_strat/(B_strat+C_strat) < 0.15)
stratParm.df$orderDiff_mh <- round(log(stratParm.df$m)-log(stratParm.df$h), 1)
ggplot(subset(stratParm.df, validPools_opt & validPools_strat)) + 
  geom_point(aes(x=cue_opt, y=cue_strat)) + 
  geom_abline(slope=1, intercept=c(0), color='grey') +
  labs(x='Optimum CUE (max B w/o compition)', y='Strategic CUE (resistant to invasion)')

#test.parm <- as.list(stratParm.df[1,])
#attach(test.parm)
#detach(test.parm)

Explore drivers of strategic and optimum differences

diff.df <- subset(stratParm.df, validPools_opt & validPools_strat, 
                  select=setdiff(names(stratParm.df), c('orderDiff_mh')))

diff.df <- melt(diff.df, measure.vars=names(parm.ls))

ggplot(diff.df) + geom_violin(aes(x=as.factor(value), y=(cue_opt-cue_strat)/cue_opt)) + 
    facet_wrap(~variable, scale='free')

We are particullary interested in differences in inputs

validParms <- subset(stratParm.df, validPools_strat & validPools_opt)
I.arr <- seq(0.1, 10, by=0.3)
Igrad <- ddply(unique(validParms[,c("b", "vmax", "kmin", "k_v_slope", "h", "m")]), 
               c("b", "vmax", "kmin", "k_v_slope", "h", "m"), function(xx){
                 ans1 <- data.frame()
                 for(I in I.arr){
                   opt_results <- optimum_cue(xx$b, xx$vmax, xx$kmin, 
                                              xx$k_v_slope, xx$h, xx$m, I)
                   names(opt_results) <- paste(names(opt_results), '_opt', sep='')
                   strat_results <- strategic_cueI(xx$b, xx$vmax, xx$kmin, 
                                                   xx$k_v_slope, xx$m, xx$h, I)
                   ans1 <- rbind.fill(ans1, 
                                      as.data.frame(c(xx[c("b", "vmax", "kmin", 
                                                           "k_v_slope", "h", "m")], 
                                                      I=I, opt_results, strat_results)))
                   
                 }
                 #cat('dim ans1:', dim(ans1), '\n')
                 ans2 <- data.frame()
                 for(I in I.arr*rnorm(length(I.arr), mean=1, sd=0.05)){
                   opt_results <- optimum_cue(xx$b, xx$vmax, xx$kmin, xx$k_v_slope, 
                                              xx$h, xx$m, I)
                   names(opt_results) <- paste(names(opt_results), '_opt', sep='')
                   strat_results <- strategic_cueI(xx$b, xx$vmax, xx$kmin, xx$k_v_slope, 
                                                   xx$m, xx$h, I)
                   ans2 <- rbind.fill(ans2, 
                                      as.data.frame(c(xx[c("b", "vmax", "kmin", 
                                                           "k_v_slope", "h", "m")], 
                                                      I=I, opt_results, strat_results)))
                 }
                 #cat('dim ans2:', dim(ans1), '\n')
                 if(nrow(ans1) != nrow(ans2)){
                   print(xx[c("b", "vmax", "kmin", "k_v_slope", "h", "m"),])
                   return(NULL)
                 }
                 
                 dCUEdI <- data.frame(dcue_opt = (ans1$cue_opt - ans2$cue_opt),
                                      dcue_strat = (ans1$cue_strat-ans2$cue_strat),
                                      dB_opt = ans1$B_opt-ans2$B_opt,
                                      dC_opt = ans1$C_opt-ans2$C_opt,
                                      dB_strat = ans1$B_strat - ans2$B_strat,
                                      dC_strat = ans1$C_strat - ans2$C_strat,
                                      dI = (ans1$I-ans2$I))
                 #cat('dim dCUEdI:', dim(dCUEdI), '\n')
                 return(cbind(ans1, dCUEdI))
               })

Igrad$validParam <- with(Igrad, B_opt+C_opt > 10 & 
                           B_opt+C_opt < 500 & 
                           B_opt/(B_opt+C_opt)> 0.001 & 
                           B_opt/(B_opt+C_opt) < 0.15 &
                           B_strat+C_strat > 10 & 
                           B_strat+C_strat < 500 & 
                           B_strat/(B_strat+C_strat)> 0.001 & 
                           B_strat/(B_strat+C_strat) < 0.15)
plot.df <- melt(Igrad, measure.vars=names(Igrad)[grepl('_(opt)|(strat)', names(Igrad))])
temp <- as.data.frame(matrix(unlist(strsplit(as.character(plot.df$variable), '_|\\.')), ncol=2, byrow=TRUE))
names(temp) <- c('variable', 'type')
plot.df$variable <- NULL
plot.df <- cbind(plot.df, temp)
plot.df$isGrad <- grepl('d', as.character(plot.df$variable))

ggplot(subset(plot.df, validParam & isGrad)) + geom_hex(aes(x=I, y=value/dI)) + facet_grid(variable~type, scales='free') + scale_y_log10()
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 70820 rows containing non-finite values (stat_binhex).

ggplot(subset(plot.df, validParam & !isGrad)) + geom_hex(aes(x=I, y=value)) + facet_grid(variable~type, scales='free')

Numerical validation code

Let’s pick some points where the gap between the Optimum CUE and Strategic CUE is small (\(<20\) percent) and large (\(>20\) percent).

bigGap.df <- subset(stratParm.df, validPools_opt & validPools_strat & abs(cue_strat-cue_opt) > 0.2)
smallGap.df <- subset(stratParm.df, validPools_opt & validPools_strat & abs(cue_strat-cue_opt) < 0.2)

numParm.df <- rbind.fill(bigGap.df[sample.int(nrow(bigGap.df), size=2),],
                         smallGap.df[sample.int(nrow(smallGap.df), size=2),])
numParm.df$index <- 1:nrow(numParm.df)
pander(numParm.df[,c(names(parm.ls), c('cue_opt','C_opt', 'B_opt', 'cue_strat','C_strat', 'B_strat'))])
Table continues below
b vmax kmin k_v_slope h m I cue_opt C_opt B_opt
-1 10 100 100 0.1 1e-04 0.1 0.9394 56.45 0.8864
1 10 100 10 1 0.001 10 0.899 314.4 8.707
-10 1 10 100 0.01 0.001 0.1 0.2828 23.72 2.158
0.1 10 100 1000 0.1 1e-04 1 0.9798 289.5 9.514
cue_strat C_strat B_strat
0.697 22.91 0.681
0.6162 71.39 6.118
0.1717 13.42 1.487
0.8687 134.8 8.57

Validate optimium CUE

opt.validate <- ddply(numParm.df, c('index', 'C_opt', 'B_opt'), function(xx){
  cue <- seq(0, 1, length=100)
  ans <- steadyState(xx$b, xx$vmax, xx$kmin, xx$k_v_slope, xx$h, xx$m, xx$I, cue)
  return(data.frame(cue=cue, C=ans$C, B=ans$B))
})

plot.df <- melt(subset(opt.validate, B>0 & C>0), measure.vars=c('C', 'B'))
ggplot(plot.df) + geom_line(aes(x=cue, y=value, group=index)) + 
  geom_vline(data=numParm.df, aes(xintercept=cue_opt)) + 
  facet_grid(variable~index, scales='free')

Validate strategic CUE

compitition.model <- function(t, y, parms){
  C <- y[1]; Bn <- y[2]; Bi <- y[3]

  ans <- with(parms,
              c(dC = I - m*C - C*(vn*Bn/(kn+C) + vi*Bi/(ki+C)),
                   dBn = cuen*vn*Bn*C/(kn+C) - hn*Bn,
                   dBi = cuei*vi*Bi*C/(ki+C) - hi*Bi)
       )
  return(list(ans))
}

runInvasions <- ddply(numParm.df, c('index'), function(xx){
  
  parm <- as.list(xx[, c('b', 'vmax', 'kmin', 'k_v_slope', 'm', 'I')])
  
  parm$hn <- xx$h
  parm$hi <- xx$h
  
  cueCombo <- expand.grid(cuen=seq(1, 100, length=50)/100,
                          rel_cuei=rnorm(10, sd=0.1))
  cueCombo$cuei <- cueCombo$cuen+rnorm(nrow(cueCombo), sd=0.1)
  cueCombo <- cueCombo[cueCombo$cuei > 0,]
  
  cueCombo$vn <- with(parm, cue_v_tradeoff(b=b, vmax=vmax, cueCombo$cuen))
  cueCombo$kn <- with(parm, v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=cueCombo$vn))
  
  cueCombo$vi <- with(parm, cue_v_tradeoff(b=b, vmax=vmax, cueCombo$cuei))
  cueCombo$ki <- with(parm, v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=cueCombo$vi))
  
  invade.df <- ddply(cueCombo, c('cuen', 'cuei'), function(cuePairs){
    comboParm <- c(parm, cuePairs)
    preInvade <- with(comboParm, steadyState(b, vmax, kmin, k_v_slope, h=hn, m, I, cue=cuen))
    if(all(preInvade > 0)){
      y0 <- list(C=preInvade$C, Bn=preInvade$B, Bi=preInvade$B*0.1)
      invasion <- lsoda(y=unlist(y0), times=c(1, 7, 30, 365, 365*10), func=compitition.model, parms=comboParm)
      names(y0) <- paste(names(y0), '0', sep='')
      if(nrow(invasion) == 5){
        finalPools <- as.data.frame(c(comboParm, y0, as.list(invasion[nrow(invasion),])))
      }else{
        finalPools <- as.data.frame(comboParm)
      }
    }else{
      finalPools <- as.data.frame(comboParm)
    }
  })
  
  return(invade.df)
})
ggplot(runInvasions) + geom_point(aes(x=cuen, y=cuei, color=Bn > Bi)) + 
  geom_vline(data=numParm.df, aes(xintercept=cue_strat)) + 
  labs(x='Native CUE', y='Invader CUE', title='Winner of 10 year invasion') +
  facet_wrap(~index)